library(SDMTools)

work.dir = "H:/Barra Work Directory/current.esoclim/"; setwd(work.dir)

pnts=cbind(x=c(112,116,116,112), y=c(-11,-11,-18.5,-18.5))
cols = colorRampPalette(c('tan', 'lightblue', 'slateblue', 'blue', 'blue4', 'black'))(101)


sum_pr =NULL
for (ii in 1:12) { cat(ii,'\n')
      tasc = read.asc.gz(paste("pr", sprintf('%02i',ii), ".asc.gz",sep=''))
      if (length(sum_pr) == 0) { #if there are no values , populate it
          sum_pr = tasc
        } else { #if already has values, add the new values to it
          sum_pr = sum_pr + tasc
		 } 
}

zlim=range(sum_pr, na.rm=T) #given what you know about zlim

png(paste('sum_pr.png',sep=''), width=7, height=7, units='cm', res=300, pointsize=5, bg='white')
image(sum_pr, zlim=zlim, ann=FALSE,axes=FALSE,col=cols)


legend.gradient(pnts,cols=cols,limits=round(zlim), title='Precipitation', cex=1)

#the new "sum_pr" layer can also be exported using the export.asc command
export.asc(x=sum_pr, file=paste(work.dir,"sum_pr_esoclim.asc",sep=""))	  

dev.off()
